多重比較(2)

Torsten Hothorn, Frank Bretz and Peter Westfall (2008), Simultaneous Inference in General Parametric Models. Biometrical Journal, 50(3), 346-363; See vignette("generalsiminf", package = "multcomp").

人工データ

set.seed(300)
count <- c(rpois(30, 10) + floor(rnorm(30, 0, 3)), rpois(30, 10) + 
    floor(rnorm(30, 0, 3)), rpois(30, 14) + floor(rnorm(30, 0, 3)))
treat <- factor(c(rep("A", 30), rep("B", 30), rep("C", 30)))
d <- data.frame(count, treat)
head(d)

##   count treat
## 1    14     A
## 2    14     A
## 3    14     A
## 4     7     A
## 5     7     A
## 6    18     A

ポアソン回帰

res <- glm(count ~ treat, d, family = poisson)
summary(res)

## 
## Call:
## glm(formula = count ~ treat, family = poisson, data = d)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -4.090  -0.923  -0.032   0.940   2.991  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.2756     0.0585   38.88   <2e-16 ***
## treatB        0.0370     0.0820    0.45     0.65    
## treatC        0.3847     0.0759    5.07    4e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 207.47  on 89  degrees of freedom
## Residual deviance: 174.69  on 87  degrees of freedom
## AIC: 555.9
## 
## Number of Fisher Scoring iterations: 4
## 

Tukeyの方法で多重比較

library(multcomp)

## Loading required package: mvtnorm

## Loading required package: survival

## Loading required package: splines

res.t <- glht(res, linfct = mcp(treat = "Tukey"))
summary(res.t)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: glm(formula = count ~ treat, family = poisson, data = d)
## 
## Linear Hypotheses:
##            Estimate Std. Error z value Pr(>|z|)    
## B - A == 0   0.0370     0.0820    0.45     0.89    
## C - A == 0   0.3847     0.0759    5.07   <1e-05 ***
## C - B == 0   0.3477     0.0750    4.63    1e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## (Adjusted p values reported -- single-step method)
## 

同時信頼区間による多重比較

res.ci <- confint(res.t)
res.ci$confint <- exp(res.ci$confint)
plot(res.ci, xlab = "relative risks")
abline(v = 1, lty = "dashed")